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ABSTRACT 

Sequential Monte Carlo algorithms (also known as particle 
filters) are popular methods to approximate filtering (and 
related) distributions of state-space models. However, they 
converge at the slow 1/y/N rate, which may be an issue 
in real-time data-intensive scenarios. We give a brief out¬ 
line of SQMC (Sequential Quasi-Monte Carlo), a variant of 
SMC based on low-discrepancy point sets proposed by jl|, 
which converges at a faster rate, and we illustrate the greater 
performance of SQMC on autonomous positioning problems. 

Index Terms — Low-discrepancy point sets; Particle fil¬ 
tering; Quasi-Monte Carlo 


show that, in several scenarios, the faster convergence of 
SQMC does more than compensate its slower running time 
and, consequently, for a given computational budget, SQMC 
typically achieves a significantly smaller error size than SMC. 

In this paper we propose to apply SQMC to the problem 
of autonomous positioning of a vehicle moving along a two 
dimensional space where, following (5|, we assume that the 
Markov transition is non Gaussian. Our numerical study show 
that for this real time application SQMC provides a much 
more accurate estimation of the position of the vehicle than 
SMC. 

2. SEQUENTIAL QUASI-MONTE CARLO 


1. INTRODUCTION 


2.1. Background on sequential Monte Carlo 


Many problems in signal processing (and related fields) can 
be formalised as the filtering of data (y t ) to recover an un¬ 
observed signal (x t ) that follows a state-space model. For 
non-linear and/or non-Gaussian state-space models, particle 
filtering | 2p| , also known as Sequential Monte Carlo (SMC), 
is now the standard approach to perform filtering; see e.g. |4j. 
However, a potential drawback of SMC for real time applica¬ 
tions is its slow 1 /y/N convergence rate (based on N simula¬ 
tions, or ‘particles’). In real time problems, the running time 
per iteration of the filtering algorithm is bounded by the time 
interval between successive observations and, consequently, 
this slow convergence rate implies that in some settings the 
approximation error of SMC might be non negligible. 

Recently, [lj proposed and studied the sequential quasi- 
Monte Carlo (SQMC) algorithm, which is a quasi-Monte 
Carlo (QMC) version of particle filtering. Based on N parti¬ 
cles, SQMC has the advantage to converge at rate o (1 / y/N ), 
i.e. at a faster rate than SMC; see Theorem 7 of |11. On the 
other hand, SQMC requires O(NlogN) operations and is 
thus slower than SMC, which has complexity O(N). But |1] 
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To introduce SMC we consider the following generic state- 
space model, described in term of probability density func¬ 
tions: 

fyt|x t ~ / y (y*|x t ), t> o 

< x t |x t _i ~ /- Y (x t |x t _i), t>i (l) 

[x 0 ~ /a Y (x o) 

where (x t )t>o is the unobservable Markov process on X C 
and (y*)t>o is the observation process. 

The typical quantity of interest in state-space models is 
the filtering distribution, that is, the distribution of x t given 
all the available observations at time t, which is given by 


p(xt|y 0:t ) = 

At 


'X t 


fo (x 0 ) n f x (x s |x a _i) / F (y s |x s )dx 0:t -1 


( 2 ) 


s=0 


where Z t is a normalising constant. Except in linear Gaus¬ 
sian models, the integrals in |2| are not tractable, but one 
may instead run a particle filter to sequentially approximate 

p(x t |yo:t). 

The basic idea of particle filtering is to use the Markov 
transition / A (x t |x t _i) to propagate the discrete approxima- 



tion (for t > 1) 

N 

P N (x t _l|y 0 : t _l) = Y 

n—1 

N 

with Y w t -1 = !, W'7-i > 0 

n—1 

of the filtering distribution at time t — 1 to the approximation 

P iV ( x t—l:f|yO:t-l) =p Ar (xt-l|y 0 :t-l)/ X (x t |x t _i) (3) 

of p(xt_i :t |yo :t _i). Then, the marginal distribution of x t 
with respect to 

P N (Xt-l:t\yO:t)cCp N (Xt-l:t\yO:t-l)f' i (yt| X t) (4) 

may be used as an approximation of the filtering distribution 
at time t. Thus, one can perform an importance sampling 
step, with proposal distribution ([3) and target distribution (|4|, 
to get a weighted particle system {VF/ 1 , x”}^ =1 which is ap¬ 
proximately distributed from p(x t |yo : t); see Algorithm [l] for 
a more precise description of particle filtering. 


Algorithm 1 SMC Algorithm (Boostrap filter) 

Operations must be performed for all n £ 1 : N 
Sample Xg from /^(x 0 ) and compute Wq = 

/ r (yo|x^)/E™=i / y ( yolx™) 

for t = 1,..., T do 

Sample a"_i from A^FF/i^), the multinomial distri¬ 
bution that produces outcome m with probability W™ 1 

Sample x” from f x (x t \xYi) and compute W" = 

/ r (y t |x?)/E™=i/ F (y t |xD 

end for 


2.2. Background on quasi-Monte Carlo 

Loosely speaking, a QMC point set u 1:Ar in [0, l) rf is a set of 
(deterministic) points which are “more uniformly” distributed 
than uniform random variates. The most classical measure of 
uniformity in the QMC literature is the so called star discrep¬ 
ancy, defined by 


which explicitly links the integration error and the equidistri- 
bution property of the point set at hand, because the quantity 
V (ip) only depends on the integrand ip; see e.g. Chap. 5 of |61 
for a definition of V(ip). 

A useful variant to QMC is randomised QMC (RQMC), 
which combines the advantages of random sampling and 
of QMC strategies. A RQMC point set u 1: v is such that 
u” ~ U([ 0,l) d ) for all n G 1 : N and D*( u 1:JV ) = 
0(N~ l (\ogN) d ) with probability one. A particularly in¬ 
teresting construction of RQMC point sets is the nested 
scrambled method for (t,m,s )-nets (see e.g. |7|, Chap. 4, 
for a definition) proposed by |8j, which allows to approxi¬ 
mate the integral of smooth functions with an error of size 
0(N~ 1 - 5+e ) for any e > 0 In addition, and contrary to 
plain QMC, no smoothness assumptions on the integrand ip 
are needed for scrambled net quadrature rules to outperform 
Monte Carlo integration [9j. This last point is particularly 
important in the context of SMC because the resampling step 
(Step [4] of Algorithm [TJ introduces discontinuities which can 
not be efficiently handled by deterministic QMC strategies. 

2.3. Sequential quasi-Monte Carlo 

The basic idea of SQMC is to replace the sampling step from 
the proposal distribution Q by a low discrepancy point set 
with respect to the same distribution. 

The classical way to transform a low discrepancy point 
set with respect to the uniform distribution (i.e. a QMC point 
set) into a low discrepancy point set with respect to a non- 
uniform distribution 7r(x) on X C is to use the inverse of 
the Rosenblatt transformation of 7r, defined by 

F^x) = (ui,... ,u d ) T , x = (xi,... ,x d ) T G X, 

where, u\ = F n i(x i), being the CDF of the marginal 
distribution of the first component (relative to tt), and for 
i > 2, Ui = F Vti (xi |xi;i_i), F W|i (-|xi :i _i) being the CDF 
of component x,;, conditional on (xi,..., Xj_i), relative to n. 

Following this idea, and assuming for the moment that 
the state variable x t is univariate, one can generate a low dis¬ 
crepancy point set (x^i, x j :Ar ) from ([3]» as follows: let u^ :Ar 
be a (R)QMC point set in [0, l) 2 , with u™ = and 

compute 


D*(u vn ) 


sup 

beQ.i)* 


i£i(u» e[ o >6 ])-nfc 

n =1 2—1 


where b = (bi ,... ,bd)- We say that u 1:JV is a QMC point set 
if D*(u 1:N ) = OiN-^logNY). 

The main motivation for using low discrepancy point sets 
in numerical integration is the Koksma-Hlawka inequality: 


N 


[ <F(u)du 
iV - J [o,i) d 


n= 1 


< V(p)D *(u 


1:JV\ 


~n _ p-1 
t_l ^P N (xt_l |yO:t' 


(«?). 


= Fjxt I -n Av?). 

i) V tJ 


However, the extension of this approach to d > 1 is not trivial 
because the distribution p Af (x t |yo : t)dx f is then a (weighted) 
sum of Dirac measures over R d . 

To overcome this difficulty, |T]| proposes to transform the 
multivariate (discrete) distribution p A, (x t _i|yo : t-i)dx t into 
a univariate (discrete) distribution p^ (ht\yo:t)dh t on [0,1) 
using the following change of variable 


x G X i->- h o ^> t (x) G [0,1], 









n = 1 n — 2 n — 3 n — 4 n — 5 



Fig. 1. The Hilbert curve is a [0,1] —► [0, l] d continuous frac¬ 
tal map, which is obtained as the limit of sequence (H n ), the 
first elements of which are represented above (for d = 2). 
Source: Marc van Dongen 


Algorithm 2 SQMC Algorithm (Boostrap filter) 

1: Operations must be performed for all n G 1 : N 

2: Generate a QMC point set Uq W in [0, l) d 

3: Compute Xg = FJ* (lift) and Wfi = 

/ v ’(yol*ft)/1!=! f Y (yo\*o) 

4: for t = 1,... ,T do 

5: Generate a QMC point set :N in [0, l) d+1 , let u” = 

«, vft), with [o,i), vfte [o,i) d 
6: Find permutation r such that < ... < 

7: Hilbert sort: find permutation at -1 such that 


where h : [0, l] d —> [0,1] is the (generalised) inverse of the 
Hilbert space filling curve H : [0,1] —> [0, l] d , and tpt ■ X —> 
[0, \} d is some user-chosen discrepancy-preserving bijection 


3.2 


between A 7 and C [0, l] d . See @ and Section 

below for more details on how to choose %j) t , and see Figure 
2.3 for a depiction of the Hilbert curve in two dimensions. 

Using this change of variable, we can see iteration t of 
SMC as an importance sampling step form 


Ph {h t -i,x t \y 0 :t-i) = 


p%(ht-^you-i)/* (xtl-ffOt-t)) (5) 


pX 


to 


Ph (h t -i,x t \y 0 :t) OC p%{ht-i,Xt\yo:t-i)f (y«|x t ) 

and we can therefore generate a low discrepancy point set 
(hi Xj :Ar ) from ([5]> as follows: let :N be a (R)QMC point 
set in [0, l) d+1 , with u" = («". v"), and compute 


ftoy) i _ 1 (x t cr l 1 1<1) ) < ... < h O 

8: Compute a]l_ 1 = where = 

Eti W:i-^ n \(n<m) 

9: Compute xf 1 = F _1 a „ (vl (n h and 

f x (- Kii 1 ) 

wp = f Y {y t W)/T,Z=if Y (yt\*?) 

10 : end for 


details on scrambling), as this makes it possible to evaluate 
the numerical error through repeated runs. 

Finally, while we presented SQMC in this specific case 
where particles are mutated according to f x (x t |x t _i), the 
Markov transition of the considered model, it of course ex¬ 
tends directly to situations where particles are mutated ac¬ 
cording to some other kernel q t (x t \x t -i) (assuming that the 
particles are reweighted accordingly, as in standard SMC). 




x?_! = H(hU), 


3. APPLICATION: AUTONOMOUS POSITIONING 


c n — U -1 (v 11 ) 


See Algorithm|2.3|for a pseudo-code description of SQMC. 


2.4. Practical implementation 


The complexity of Algorithm |2.3| is 0(N log iV), because 
it performs two sorting steps at each iteration. Regarding 
the practical implementation of Algorithm [23] note that: (a) 
QMC generation (Steps 2 and 5) routines are available in 
most software (e.g. package randtoolbox in R, or class 
qrandset in the Statistics toolbox of Matlab); to com¬ 
pute the a™_ 1 ’s (Step 8), one may use the standard approach 
based on sorted uniforms for resampling; and (c) in order to 


compute h, see e.g. 1101, and Chris Hamilton’s C++ program 
available at https://web.cs.dal.ca/~chamilto/ 
hilbert/index.html. 

Our SQMC implementation is available at https:// 
bitbucket. org/mgerber/sqmc We shall use RQMC 
(randomised QMC) point sets in our simulations (more pre¬ 
cisely scrambled Sobol’ sequences; see |9 11 121 for more 


3.1. Model description 

We consider the problem of autonomous positioning of a ve¬ 
hicle moving in a two dimensional space. To determine its 
position, the vehicle estimates its speed every T s > 0 seconds 
and measures the power of dy > 1 radio signals. We sup¬ 
pose that the radio signals are emitted from known locations 
rt G M 2 , i = 1,... ,d y , and that the corresponding attenua¬ 
tion coefficients at are known as well. This positioning prob¬ 
lem admits the following state space representation (see G3 
and (5)) 

f y ti = 10log 10 ( ii^xjc, ) + vu, t> 0 

< xt = x t _i + T s v t + T s e t , t> 1 ( 6 ) 

[x 0 ~ A4(0, 1 2 ) 


where * G 1 : d y , x t £ R 2 is the position of the vehicle at 
time t, v t is a measure of its speed, which is assumed to be 
constant over successive time intervals of T s seconds, e t and 
i/ y = (vu, ■ ■ ■, t) represent measurement errors while y it 
is the power received at time t by emitter i. In the sequel. Pot 



























and thus, a reasonable choice for x ti and xu is 

t 

X ti , Xu = ^2v s ± 2 v /Var(x 0i ) + t T s 2 Var(e H ). 

8 =o 

Simulation results are presented for N £ {2 s ,..., 2 16 }, 
where 2 is the base of the Sobol’ sequence. Taking a power 
of 2 for the number of simulations is the standard approach 
in QMC integration based on Sobol’ sequence because both 
good theoretical and empirical results are obtained for this 
choice of N. However, this restriction is non necessary for 
QMC to outperform Monte Carlo methods and little gain may 
be expected in the context of SQMC, see 03 for more details 
on this point. 

-15000 -10000 -5000 6 



Fig. 2. Trajectory of a vehicle evolving for 15 minutes and 
starting at a location close to (0, 0). The dots show the loca¬ 
tions of the 5 emitters. 


is the initial signal from emitter i and, following |[5j, we sup¬ 
pose that all the error terms are independent and distributed 
according to a Laplace distribution with parameter 0.5. 


3.2. Simulation set-up 


To compare the performance of SQMC and SMC for this 
tracking problem we simulate the trajectory of a vehicle 
evolving during 15 minutes according to |6]). We assume that 
the sample period is T s = 1 second, that d y = 5 (5 emitters) 
and that a, = 0.95 for all i = 1,..., d y . The resulting tra¬ 
jectory and the locations of the emitters are shown in Figure 

m 


The SMC algorithm is implemented using systematic re¬ 
sampling [141, which is usually recognised as being the most 
efficient resampling strategy. 

SQMC is implemented using nested scrambled Sobol’ 
sequences for the point sets uJ :W . As described above, 
we need to use a mapping ip t to map the particles gen¬ 
erated at iteration t of SQMC into the unit square before 
performing the Hilbert sort. Following [jTJ, we chose for 
il>t a component-wise (re-scaled) logistic transform; that is, 
V>t(x) = (V’tiOn), ipn{xi)) with 


1pti{Xi) 


1 + exp 



Xi - x u \ 
Xu - X ti ) 


-1 


i = 1,2. 


and where the time varying constants x tl and x ti are used to 
solve numerical problems due to high values of \xi\. More 
precisely, these constants should be chosen such that, with 
high probability, x t i £ [x ti , x t i]- To this aims, note that 

Var(a; ti ) = Var(xoi) + t T s 2 Var(e H ). 


3.3. Results 

In Figure [3] we compare the mean square error (MSE) of the 
filtering expectation estimate obtained from SQMC and SMC, 
as a function of t, for N £ {2 8 ,2 10 , 2 16 }. To save space, 
only the results for the first component of x t are presented; 
the results for the second component are essentially the same. 
One observes that the performance gain of SQMC (relative to 
standard particle filtering) increases quickly with N. 

We now study the amount of CPU time required to 
have a “reasonable” Monte Carlo error using both SMC and 
SQMC. Letting x t be an estimate of the filtering expectation 
E[x f |y 0: t], we consider the Monte Carlo error to be reason¬ 
able if it is small compared to the posterior variance, that is, 
if MSE(ijt) < 6 2 Var(;cit|yo:t) for i = 1, 2 and where we set 

5 = 0.01. 

Figure[4]shows the number of time steps t £ {0,..., 899} 
for which this criterion is not met, as a function of the CPU 
budget (i.e. CPU time per iteration). To increase the CPU 
budget, we simply increase N. We observe that much better 
results are achieved using SQMC. Indeed, when the CPU bud¬ 
get is 0.05s per iteration, the SMC error is too large for more 
than 600 time steps, while a CPU budget of 0.07s is enough to 
estimate both coordinates of x* for all iterations with SQMC. 

4. CONCLUSION 

In this paper we have illustrated the potential of sequential 
quasi-Monte Carlo for real time signal processing process¬ 
ing problems with a non-linear and non-Gaussian state-space 
model for autonomous positioning. Compared to Monte 
Carlo particle filtering, dramatic variance reductions are ob¬ 
served when SQMC is used, both as a function of the number 
of particles and of CPU time. In real time application, the 
running time of the filtering algorithm is a crucial element 
and, concerning this point, we believe that significant im¬ 
provement can be achieved for SQMC, notably concerning 
the Hilbert sort step. For instance, the computations of the 
Hilbert indices involve only bits operations and therefore 
GPU computing may allow for dramatic cost reductions. 
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Fig. 3. Filtering of the state-space model The plot gives 
the gain factor, defined as the MSE(SMC) over MSE(SQMC), 
as a function of t, for the estimation of E[a;it|yo:t)]- The re¬ 
sults are obtained from 100 independent runs of SMC and 
SQMC. 



CPU time in second 

Fig. 4. Filtering of the state-space model & The plot 
gives the number of time steps t £ {0,..., 899} such that 
MSE(£’ t i) > 0.01 2 Var(a;ti|yo:t) as a function of the CPU 
budget (average CPU time per iteration), where Xu is either 
the SQMC (solid lines) or the SMC (dashed lines) estimate 
of E[xtj|yo;t], i = 1,2. The results are obtained from 100 
independent runs of SMC and SQMC. 
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